Method for characterizing the fracture network of a fractured reservoir and method for developing it

ABSTRACT

The invention is a method for constructing a representation of a fluid reservoir traversed by a fracture network and by at least one well. The reservoir is discretized into a set of grid cells and the fractures are characterized by statistical parameters from observations of the reservoir. An equivalent permeability tensor and an average fracture opening is constructed from an image representative of the fracture network delimiting porous blocks and fractures is then deduced from the statistical parameters. A first elliptical boundary zone centered on the well and at least a second elliptical boundary zone centered on the well which form an elliptical ring with the elliptical boundary of the first zone are defined around the well. The image representative of the fracture network is simplified in a different manner for each of the zones which is used to construct the representation of the fluid reservoir.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to development of underground reservoirs,such as hydrocarbon reservoirs comprising a fracture network. Inparticular, the invention relates to a method for characterizing thefracture network and constructing a representation of the reservoir. Therepresentation is used to optimize the management of development througha prediction of the fluid flows likely to occur through the medium tosimulate hydrocarbon production according to various productionscenarios.

2. Description of the Prior Art

The petroleum industry, and more precisely exploration and developmentof reservoirs, notably petroleum reservoirs, require knowledge of theunderground geology which is as perfect as possible so as to efficientlyprovide evaluation of reserves, production modelling or developmentmanagement. In fact, determining the location of a production well or ofan injection well, the drilling mud composition, the completioncharacteristics, selection of a hydrocarbon recovery method (such aswaterflooding for example) and of the parameters required forimplementing this method (such as injection pressure, production flowrate, etc.) requires good knowledge of the reservoir. Reservoirknowledge notably means knowledge of the petrophysical properties of thesubsoil at any point in space.

The petroleum industry has therefore combined for a long time field(in-situ) measurements with experimental modelling (performed in thelaboratory) and/or numerical modelling (using softwares). Petroleumreservoir modelling is a technical stage that is essential for anyreservoir exploration or development procedure. The goal of modelling isto provide a description of the reservoir.

Fractured reservoirs are an extreme type of heterogeneous reservoirscomprising two very different media, a matrix medium containing themajor part of the oil in place and having a low permeability, and afractured medium representing less than 1% of the oil in place andhighly conductive. The fractured medium itself can be complex, withdifferent sets of fractures characterized by their respective density,length, orientation, inclination and opening.

Those in charge of the development of fractured reservoirs need toperfectly know the role of fractures. What is referred to as a“fracture” is a plane discontinuity of very small thickness in relationto the extent thereof, representing a rupture plane of a rock of thereservoir. On the one hand, knowledge of the distribution and of thebehavior of these fractures allows optimizing the location and thespacing between wells to be drilled through the oil-bearing reservoir.On the other hand, the geometry of the fracture network conditions thefluid displacement, on the reservoir scale as well as the local scalewhere it determines elementary matrix blocks in which the oil istrapped. Knowing the distribution of the fractures is therefore alsovery helpful at a later stage to the reservoir engineer who wants tocalibrate the models which have been constructed to simulate thereservoirs in order to reproduce or to predict the past or futureproduction curves. Geosciences specialists therefore havethree-dimensional images of reservoirs allowing locating a large numberof fractures.

Thus, in order to reproduce or to predict (that is “simulate”) theproduction of hydrocarbons when starting production of a reservoiraccording to a given production scenario (characterized by the positionof the wells, the recovery method, etc.), reservoir engineers use acomputing software referred to as reservoir simulator (or flowsimulator) that calculates the flows and the evolution of the pressureswithin the reservoir represented by the reservoir model. The results ofthese computations enable prediction and optimization of the reservoirin terms of flow rate and/or of amount of hydrocarbons recovered.Calculation of the reservoir behavior according to a given productionscenario constitutes a “reservoir simulation”.

There is a well-known method for optimizing the development of a fluidreservoir traversed by a fracture network, wherein fluid flows throughthe reservoir are simulated through simplified but realistic modellingof the reservoir. This simplified representation is referred to as“double-medium approach”, described by Warren J. E. et al. in “TheBehavior of Naturally Fractured Reservoirs”, SPE Journal (September1963), 245-255. This technique considers the fractured medium as twocontinua exchanging fluids with one another: matrix blocks and fractureswhich is referred to as a “double medium” or “double porosity” model.Thus, “double-medium” modelling of a fractured reservoir discretizes thereservoir into two superposed sets of cells (referred to as grids)making up the “fracture” grid and the “matrix” grid. Each elementaryvolume of the fractured reservoir is thus conceptually represented bytwo cells, a “fracture” cell and a “matrix” cell, coupled to one another(i.e. exchanging fluids). In the reality of the fractured field, thesetwo cells represent all of the matrix blocks delimited by fracturespresent at this point of the reservoir. In fact, in most cases, thecells have hectometric lateral dimensions (commonly 100 or 200 m)considering the size of the fields and the limited possibilities ofsimulation softwares in terms of computing capacity and time. The resultthereof is that, for most fractured fields, the fractured reservoirelementary volume (cell) comprises innumerable fractures forming acomplex network that delimits multiple matrix blocks of variabledimensions and shapes according to the geological context. Eachconstituent real block exchanges fluids with the surrounding fracturesat a rate (flow rate) that is specific thereto because it depends on thedimensions and on the shape of this particular block.

In the face of such a geometrical complexity of the real medium, theapproach is for each reservoir elementary volume (cell), in representingthe real fractured medium as a set of matrix blocks that are allidentical, parallelepipedic, delimited by an orthogonal and regularnetwork of fractures oriented in the principal directions of flow: Foreach cell, the so-called “equivalent” permeabilities of this fracturenetwork are thus determined and a matrix block referred to as“representative” (of the real (geological) distribution of the blocks),single and of parallelepipedic shape, is defined. It is then possible toformulate and to calculate the matrix-fracture exchange flows for this“representative” block and to multiply the result by the number of suchblocks in the elementary volume (cell) to obtain the flow on the scaleof this cell.

It can however be noted that calculation of the equivalentpermeabilities requires knowledge of the flow properties (that is theconductivities) of the discrete fractures of the geological model.

It is therefore necessary, prior to constructing this equivalentreservoir model (referred to as “double-medium reservoir model”) asdescribed above, to simulate the flow responses of some wells (transientor pseudo-permanent flow tests, interferences, flow measurement, etc.)on models extracted from the geological model giving a discrete(realistic) representation of the fractures supplying these wells.Adjustment of the simulated pressure/flow rate responses on the fieldmeasurements allows the conductivities of the fracture families to becalibrated. Although it covers a limited area (drainage area) around thewell only, such a well test simulation model still comprises a verylarge number of calculation nodes if the fracture network is dense.Consequently, the size of the systems to be solved and/or thecomputation time often remain prohibitive.

SUMMARY OF THE INVENTION

To overcome this difficulty, the invention comprises a simplification ofthe fracture networks on the local scale of the well drainage area, soas to be able to simulate the fractured reservoir well tests and thus tocalibrate the conductivities of the fracture families. This hydrauliccalibration of the fractures leads to a set of parameters characterizingthe fracture network (or fracture model). This fracture model isthereafter used to construct a double-medium flow model on the reservoirscale.

The invention thus relates to a method for optimizing the development ofa fluid reservoir traversed by a fracture network and by at least onewell, wherein a representation of the fluid reservoir is constructed.The reservoir is discretized into a set of cells and the fractures arecharacterized by statistical parameters (PSF) from observations of thereservoir. The method comprises the following stages:

-   -   a) deducing from the statistical parameters (PSF) an equivalent        permeability tensor and an average opening for the fractures,        from which an image representative of the fracture network        delimiting porous blocks and fractures is constructed;    -   b) deducing from the tensor a direction of flow of the fluid        around the well;    -   c) defining around the well a first elliptical boundary zone        centered on the well and containing the well, and at least a        second elliptical boundary zone centered on the well and forming        an elliptical ring with the elliptical boundary of the first        zone, the zones being oriented in the direction of flow of the        fluid;    -   d) simplifying the image representative of the fracture network        in a different manner in each one of the zones;    -   e) using the simplified image to construct the representation of        the fluid reservoir; and    -   f) using the representation of the fluid reservoir and a flow        simulator to optimize the development of the fluid reservoir.

According to the invention, the statistical parameters (PSF) can beselected from among the following parameters: fracture density, fracturelength, fracture orientation, fracture inclination, fracture opening andfracture distribution within the reservoir.

According to an embodiment, an aspect ratio is determined for each zone,defined from the lengths of the axes of the ellipse making up theboundary of the zone, so as to reproduce a flow anisotropy around thewell, and the zones are constructed so as to respect the aspect ratio.This aspect ratio can be determined by the principal values of thepermeability tensor.

According to an embodiment, a distance is defined between the boundariesbetween zones, so as to give an equal weight to each zone in terms ofpressure difference recorded in each zone under permanent flow regimeconditions. This distance can be defined by setting the lengths of oneof the two axes of two successive ellipsis at values in geometricprogression of constant ratio.

According to another embodiment, three zones are constructed, a firstzone (ZNS) containing the well, wherein no image simplification isprovided, a second zone (ZP) in contact with the first zone, wherein afirst image simplification is carried out, and a third zone (ZL) incontact with the second zone, wherein a second image simplification iscarried out with the second simplification being more significant thanthe first one.

Advantageously, the second and third zones can be divided into sub-zonesby applying the following stages:

the second zone is divided into a number of sub-zones equal to a numberof blocks of cells present in the zone with a block of cells designatinga vertical pile of cells; the third zone is divided by carrying out thefollowing stages:

-   -   dividing every degree of the boundary of the third zone into 360        arcs;    -   defining a sub-zone by connecting end points of each of the arcs        to the center of the ellipse forming the boundary;    -   for each one of the sub-zones, calculating an equivalent        fracture permeability tensor from which an orientation of the        flows in the sub-zone is determined;    -   comparing the equivalent fracture permeability values and the        flow orientation between neighboring sub-zones; and    -   grouping neighboring sub-zones together into a single sub-zone        when a difference between the permeability values is below a        first threshold and when a difference between the flow        orientations is below a second threshold.

According to the invention, the image can be simplified by carrying outthe following stages:

-   -   constructing a fracture network equivalent of the image (RFE),        by means of a so-called Warren and Root representation, wherein        the network is characterized by fracture spacings (s₁ ^(fin), s₂        ^(fin)) in two orthogonal directions of principal permeability,        by a fracture opening parameter (e^(fin)), by fracture        conductivities (C_(f1) ^(fin) and C_(f2) ^(fin)) and a        permeability (k_(m) ^(fin)) of a matrix medium between        fractures; and    -   simplifying the equivalent fracture network (RFE) by a network        fracture spacing coefficient (G) whose value is less than a        value G_(max-zone) defined on each of the zones in order to        guarantee sufficient connectivity between simplified zones and        non-simplified zones.

Value G_(max-zone) can be defined as follows for a given sub-zone:

$G_{\max - {zone}} = \frac{DLM}{6 \cdot {{Max}\left( {s_{1}^{fin},s_{2}^{fin}} \right)}}$with:

-   -   DLM: minimum lateral dimension of the given sub-zone;    -   s₁ ^(fin), s₂ ^(fin): fracture spacings in the so-called Warren        and Root representation.

Finally, the invention also relates to a method for optimizing themanagement of a reservoir. It comprises the following stages:

-   -   repeating stage a) while modifying the statistical parameters        (PSF) so as to minimize a difference between a well test result        and a well test simulation result from the simplified image;    -   associating with each one of the cells at least one equivalent        permeability value and an average opening value for the        fractures with the values being determined from the modified        statistical parameters (PSF);    -   simulating fluid flows in the reservoir with a flow simulator,        and the equivalent permeability values and the average opening        values of the fractures associated with each one of the cells;    -   selecting a production scenario allowing optimizing the        reservoir production with the fluid flow simulation; and    -   developing the reservoir according to the scenario allowing        optimizing the reservoir production.

BRIEF DESCRIPTION OF THE DRAWINGS

Other features and advantages of the method according to the inventionwill be clear from reading the description hereafter of embodimentsgiven by way of non limitative example, with reference to theaccompanying figures wherein:

FIG. 1 illustrates the various stages of the method according to theinvention;

FIG. 2 illustrates a realization of a fracture/fault network on thereservoir scale;

FIG. 3 illustrates an initial discrete fracture network (DFN);

FIG. 4 illustrates a so-called Warren and Root equivalent fracturenetwork (RFE);

FIG. 5 illustrates the creation of zones and sub-zones necessary forsimplification of the equivalent fracture network (RFE);

FIG. 6 illustrates a simplified equivalent fracture network (RFES)according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

The method according to the invention for optimizing the development ofa reservoir using the method of the invention for characterizing thefracture network comprises four stages, as illustrated in FIG. 1:

1—Discretization of the reservoir into a set of cells (MR)

2—Modelling of the fracture network (DFN, RFE, RFES)

3—Simulation of the fluid flows (SIM) and optimization of the reservoirproduction conditions (OPT)

4—Optimized (global) development of the reservoir (EXPLO)

1—Discretization of the Reservoir into a Set of Cells (MR)

The petroleum industry has combined for a long time field (in-situ)measurements with experimental modelling (performed in the laboratory)and/or numerical modelling (using softwares). Petroleum reservoirmodelling thus is an essential technical stage with a view to reservoirexploration or development. The goal of such modelling is to provide adescription of the reservoir, characterized by the structure/geometryand the petrophysical properties of the geological deposits orformations therein.

These modellings are based on a representation of the reservoir as a setof cells. Each cell represents a given volume of the reservoir and makesup an elementary volume of the reservoir. The cells in their entiretymake up a discrete representation of the reservoir which is referred toas geological model.

Many software tools are known allowing construction of such reservoirmodels from data (DG) and measurements (MG) relative to the reservoir.

FIG. 2 illustrates a two-dimensional view of a reservoir model. Thefractures are represented by lines. The cells are not shown.

2—Modelling the Fracture Network

In order to take into account the role of the fracture network in thesimulation of flows within the reservoir, it is necessary to associatewith each of these elementary volumes (cells of the reservoir model) amodelling of the fractures.

Thus, one object of the invention relates to a method for constructing arepresentation of a fluid reservoir traversed by a fracture network andby at least one well. This method comprises discretizing the reservoirinto a set of cells (stage 1 described above). The method then comprisesthe following stages:

-   a. characterizing the fractures by statistical parameters (PSF) from    reservoir observations;-   b. determining from these statistical parameters (PSF) an equivalent    permeability tensor and an average opening of the fractures, from    which an image representative of the fracture network delimiting    porous blocks and fractures is constructed;-   c. determining from the tensor a direction of flow of the fluid    around the well;-   d. defining around the well a first elliptical boundary zone    centered on the well and containing the well and at least a second    elliptical boundary zone centered on the well within an inner    boundary merging with the elliptical boundary of the first zone with    the zones being oriented in the direction of flow of the fluid;-   e. simplifying the image representative of the fracture network in    each cell belonging to at least one zone;-   f. repeating b) while modifying the statistical parameters (PSF) to    minimize the difference between the well test result and the well    test simulation result from the simplified image; and-   g. associating with each cell at least one equivalent permeability    value and a fracture average opening value with these values being    determined from the modified statistical parameters (PSF).

These stages are described in detail hereafter.

Fracture Characterization

The statistical reservoir characterization is based upon carrying outdirect and indirect reservoir observations (OF). This characterizationuses 1) well cores extracted from the reservoir, on which a statisticalstudy of the intersected fractures is performed, 2) outcrops which arecharacteristic of the reservoir which has the advantage of providing alarge-scale view of the fracture network, and 3) seismic images allowingthe identification of major geological events.

These measurements allow characterizing the fractures by statisticalparameters (PSF) which are their respective density, length,orientation, inclination and opening, and their distribution within thereservoir.

At the end of this fracture characterization stage, statisticalparameters (PSF) are obtained describing the fracture networks fromwhich realistic images of the real (geological) networks can bereconstructed (generated) on the scale of each cell of the reservoirmodel considered (simulation domain).

The goal of characterization and modelling of the reservoir fracturenetwork is to provide a fracture model validated on the local flowsaround the wells. This fracture model is then extended to the reservoirscale in order to achieve production simulations. Flow properties aretherefore associated with each cell of the reservoir model (MR)(permeability tensor, porosity) of the two media (fracture and matrix).

These properties can be determined either directly from the statisticalparameters (PSF) describing the fracture networks, or from a discretefracture network (DFN) obtained from the statistical parameters (PSF).

Constructing a Discrete Fracture Network (DFN)·FIGS. 2 and 3

Starting from a reservoir model of the reservoir being studied, adetailed representation (DFN) of the internal complexity of the fracturenetwork which is as accurate as possible in relation to the direct andindirect reservoir observations, is associated with each cell. FIG. 2illustrates a realization of a fracture/fault network on the scale of areservoir. Each cell of the reservoir model thus represents a discretefracture network delimiting a set of porous matrix blocks, of irregularshape and size, delimited by fractures. Such an image is shown in FIG.3. This discrete fracture network constitutes an image representative ofthe real fracture network delimiting the matrix blocks.

Constructing a discrete fracture network in each cell of a reservoirmodel can be achieved using known modelling softwares, such as theFRACAFlow® software (IFP, France). These softwares use the statisticalparameters determined in the fracture characterization stage.

The next stage determines the flow properties of the initial fractures(C_(f), e) and then calibrates these properties by well test simulationson discrete local flow models obtained from the realistic image of thereal (geological) fracture network on the reservoir scale. Although itcovers a limited area (drainage area) around the well only, such a welltest simulation model still comprises a very large number of calculationnodes if the fracture network is dense. Consequently, the size of thesystems to be solved and/or the computation time often remainprohibitive, hence the necessity to use a fracture networksimplification procedure.

Fracture Network Simplification—FIG. 5

Because of its extreme geometrical complexity, the fracture networkobtained in the previous stage, representative of the real fracturedreservoir, cannot be used to simulate, that is reproduce and/or predict,the local flows around the well.

In order to overcome this obstacle, the method according to theinvention uses a procedure based on the division of the simulationdomain (that is the reservoir model) into at least three types of zonearound each well (FIG. 5):

-   -   A first zone wherein no fracture network simplification is        performed. This zone contains the well and its center. It is        denoted by ZNS (non-simplified zone);    -   A second zone which is in contact with the first zone, wherein a        moderate simplification of the fracture network is performed.        This zone is denoted by ZP (zone to be simplified close to the        well);    -   A third zone, in contact with the second zone, wherein a        significant simplification of the fracture network is performed.        This zone is denoted by ZL (zone at a distance from the well).

The invention is not limited to the definition of three zones. It isalso possible to divide the domain into n zones in which thesimplification of the network increases from zone 1 (ZNS) to zone n (thefurthest from the well). It is thus possible to create a zone ZNS, n1zones of type ZP and n2 zones of type ZL.

Constructing the Zones—FIG. 5

The goal of fracture network modelling is to simulate the well flowresponses (transient or pseudo-permanent flow tests, interferences, flowmeasurement, etc.). It consists for example in simulating the oilproduction via each well drilled through the reservoir.

For each well, each zone is defined according to an exterior boundaryforming an ellipse centered on the well. The three zones are thusconcentric and with elliptical boundaries. The two simplified zones ZPand ZL have inner boundaries corresponding to the outer boundary of therespective zones ZNS and ZP. Except for the non-simplified zone, thezones thus are elliptical rings centered on the well. To construct eachzone, it is a must to define:

-   -   the orientation of the ellipse, that is the direction of the        major axis of the ellipse (perpendicular to the direction of the        minor axis); and    -   the dimensions of the ellipse, that is the length of the axes.

Their orientation is determined by the directions of flow determinedfrom a calculation of equivalent permeabilities on zone ZNS. This typeof equivalent permeability calculation is well known. It is possible touse, for example, the numerical method of calculating fractured mediaequivalent properties, implemented in the FracaFlow software (IFPEnergies nouvelles, France) and described hereafter.

According to this method, a permeability tensor representative of theflow properties of the discretized fracture network (DFN) can beobtained via two upscaling methods.

-   -   The first method which is an analytical method referred to as        local analytical upscaling, is based on an analytical approach        described in the following documents:

-   Oda M. (1985): Permeability Tensor for Discontinuous Rock Masses,    Geotechnique Vol 35, 483-495; and

-   Patent application EP 2 037 080.

It affords the advantage of being very fast. Its range of application ishowever limited to well connected fracture networks. In the oppositecase, major errors on the permeability tensor can be observed.

-   -   The second method which is a numerical method referred to as        local numerical upscaling, is described in the following        documents:

-   Bourbiaux, B., et al., 1998, “A Rapid and Efficient Methodology to    Convert Fractured Reservoir Images into a Dual-Porosity Model”, Oil    & Gas Science and Technology, Vol. 53, No. 6, November-December    1998, 785-799.

-   French Patent 2.757.947, corresponding to U.S. Pat. No. 6,023,656    for equivalent permeabilities, and French Patent 2,757,957,    corresponding to U.S. Pat. No. 6,064,944, for equivalent block    dimensions.

It is based on the numerical solution of the equations of flow on adiscrete grid of the fracture network for various boundary conditions ofthe computation block considered. The equivalent permeability tensor isobtained by identification of the ratios between flow rate and pressuredrop at the boundaries of the computation block. This approach, which ismore expensive than the previous one, has the advantage ofcharacterizing a given network (even weakly connected well).

According to an embodiment, it is possible to select one or the other ofthe two previous methods to optimize the accuracy and speed of thecomputations, by applying the method described in EP Patent Application2,037,080, based on the computation of a connectivity index.

This technique allows, in a preliminary stage, determination of thepermeability tensors of some cells of the reservoir model surroundingthe well, and is considered to be representative of the flow of ZNS.Diagonalization of these permeability tensors provides the eigenvectorsoriented in the principal directions of flow that are being sought. Itis then possible to orient the elliptical domain ZNS centered on thewell along the semi-major axis of the permeability ellipse determined bythe preliminary computations.

Then, in order to define the dimensions of the ellipse, an aspect ratiois defined for this ellipse, as well as a distance between concentricellipsis (distance between the inner and outer elliptical boundaries ofa given concentric ring):

-   -   aspect ratio of the ellipsis. By denoting by L_(max) the        half-length of the major axis of the ellipse, and by L_(min) the        half-length of the minor axis of the ellipse, this aspect ratio        is defined by L_(max)/L_(min). It should not be fixed        arbitrarily but, and to the contrary, it should be selected to        reproduce the flow anisotropy. The calculated equivalent        permeability tensor can be used to determine the orientation. If        the principal values of this tensor are denoted by K_(min) and        K_(max), then the ellipsis are oriented in the principal        permeability directions with an aspect ratio L_(max)/L_(min)        equal to the square root of ratio K_(max)/K_(min):

$\frac{L_{\max}}{L_{\min}} = \sqrt{\frac{K_{\max}}{K_{\min}}}$

-   -   which is a distance between concentric ellipsis (boundaries        between zones).        A dimensioning criterion in accordance with a modelling accuracy        uniformly distributed over the simulation domain sets the        half-lengths of the major axis L_(max) (or of the minor axis        L_(min)) of 2 successive ellipsis i+1 and i at values in        geometric progression of constant ratio r (equal to 2 for        example), that is:

${\frac{L_{\max{({i + 1})}}}{L_{\max{(i)}}} = {\frac{L_{\min{({i + 1})}}}{L_{\min{(i)}}} = r}},$for any I, with initialization to the value L_(max0) of thenon-simplified zone. This rule allows giving an equal weight to eachzone I (i=1 to n) in terms of pressure difference observed in each ringunder permanent flow regime conditions.

Once this global dimensioning is achieved, the sub-zone delimitation andsimplification procedures are implemented which are also based on themethods of calculating the equivalent properties of fractured media.

Constructing Sub-Zones in Each Zone—FIG. 5

Each zone, except zone ZNS, is then subdivided into sub-zones (ssZP,ssZL1, ssZL2) within which simplification of the fracture network isperformed. The number of sub-zones per zone depends on the type of thezone (ZNS, ZP or ZL) and on the heterogeneity thereof. Thus:

-   -   Zone ZNS is to be kept intact (non simplified); no sub-zone is        created therein,    -   The zones to be simplified that are the closest to a well (ZP        type zones) require particular attention. In order to correctly        model the local variations of the flow properties in this zone,        zones ZP are divided into a number of sub-zones equal to the        number of blocks of cells present in zones ZP. The term “block        of cells” is used to designate a “stack” of vertical cells (of        CPG type (Corner Point Grid) for example) of the reservoir        model, delimited by the same sub-vertical upright poles. They        are therefore no equivalent blocks;    -   The zones to be simplified that are the furthest from a well (ZL        type zones) are concentric elliptical rings that cover        increasingly large surfaces as the distance from the well        increases. Being further away from the wells, it is acceptable        to be less accurate regarding heterogeneities detection than in        the case of ZP type zones. The heterogeneity of the flow        properties of a ZL type zone remains however the principal        factor controlling the subdivision into sub-zones. To sample a        ZL type zone, an angular scan centered on the well is performed        on the outer elliptical limit separating the zone ZL being        considered from its outer neighbor. For every degree, a block of        cells is selected. Thus, zone ZL is characterized by at most 360        blocks. For each block, the equivalent fracture permeability        tensor is calculated. This tensor allows characterizing the        dynamic properties of the fracture network of the block being        studied. The permeability values and the flow orientation        obtained are then compared among neighboring blocks. If the        properties of two neighboring blocks are close, these two blocks        are considered to belong to the same sub-zone. In the opposite        case (20% difference on the principal permeability values or 10        degree difference on the principal permeability directions for        example), the two blocks are assigned to different sub-zones.        The outer limit of the zone ZL is considered is divided into        arcs of elliptical shape (FIG. 1). The sub-zones are then        obtained by connecting the end points of each arc to the center        (well) of the ellipse (dotted lines in FIG. 5). Each sub-zone is        thus defined as the area contained between the intra-zone radial        partition (dotted lines) and the inter-zone elliptical limits        (full line).

As mentioned above, the delimitation of the sub-zones (limit points ofthe inter-zone elliptical boundary arcs) is based on the comparison ofthe equivalent permeabilities calculated on the “neighboring” blocksmarking these arcs. The analytical method of computing the equivalentpermeabilities is preferably used because the goal is, in this case, todetermine whether a block has the same dynamic behavior as theneighboring block. Although, for a weakly connected network, theanalytical approach provides erroneous results, errors are systematic,which are similar from one block to the next, which allows comparison ofthe results between blocks which clearly does not require a highaccuracy considering the simple zone definition objective. Thus, theanalytical approach is totally justified, with the considerableadvantage of enabling must faster calculations, thus guaranteeingpractical feasibility.

Simplifying the Fracture Network in the Sub-Zones (RFE, RFES)

Once the zones ARE divided into sub-zones, the upscaling calculationsallow replacement of the fracture network of these sub-zones by asimplified network having the same flow properties as the initialnetwork. In this case, and unlike what is written above, the calculationof the equivalent fracture permeability tensor has to be as accurate aspossible.

In order to make the most of the advantages afforded by the twoupscaling methods, the permeability tensor is determined by one or theother of these two methods for example, according to the selectionprocedure described in document EP Patent 2,037,080, based on the valueof the connectivity index of the fracture network. This index, which isrepresentative of the ratio between the number of intersections betweenfractures and the number of fractures, is calculated for each unit ofthe block being considered (that is 2D). Its value allows consideringthe network as very well connected, weakly/badly connected ornon-connected. The upscaling method is then selected as follows Delorme,M., Atfeh, B., Alken, V. and Bourbiaux, B. 2008, Upscaling Improvementfor Heterogeneous Fractured Reservoir Using a GeostatisticalConnectivity Index, edited in Geostatistics 2008, VIII InternationalGeostatistics Congress, Santiago, Chile:

-   -   A well-connected network in this case is characterized by a        connectivity index close to or exceeding 3 (at least 3        intersections per fracture of the network on average). The        analytical upscaling method is selected because its accuracy is        guaranteed considering the good connectivity of the network with        the essential additional advantage of fastness;    -   A weakly/badly-connected network in this case is when the        connectivity index ranges between 1 and 3 (which corresponds to        a number of fracture intersections ranging between one and three        times the number of fractures), and the numerical upscaling        method is selected to reliably calculate the permeability        tensor;    -   When the network is very poorly or even not connected, which        occurs when number of intersections is close to or lower than        the number of fractures). The original fractures (which are few)        are kept. That is, the sub-zone in question is not simplified.

Once these equivalent permeability calculations are carried out for eachsub-zone, the other equivalent flow parameters characterizing thesimplified sub-zones are readily determined according to the followingmethods and equivalences:

-   -   Calculation of a first equivalent network (so-called Warren and        Root parallelepipedic network), which is referred to as “fine”,        directly resulting from the method of French Patent 2,757,957        corresponding to U.S. Pat. No. 6,064,944. This network (FIG. 4)        is characterized by spacings of fractures s₁ ^(fin), s₂ ^(fin)        in the 2 orthogonal directions of principal permeability 1 and 2        defining the fracture network.

An additional parameter, the opening of the fractures (e^(fin)),characterizes the fractures of this fine network. The value of thefracture openings is in practice nearly always negligible in relation tothe fracture spacing. This hypothesis is taken into account in thefollowing formulas, where the same opening value is assumed for the 2fracture families. Considering the equality of porosities of the initialnetwork DFN ( _(f)) and of the fine equivalent network, e^(fin) isdeduced from the fracture volume of the initial network DFN, V_(f)^(init), and from the total rock volume V_(T) as follows:

$\begin{matrix}{e^{fin} = {{\frac{1}{\left( {\frac{1}{s_{1}^{fin}} + \frac{1}{s_{2}^{fin}}} \right)}\frac{V_{f}^{init}}{V_{T}}} = {\frac{1}{\left( {\frac{1}{s_{1}^{fin}} + \frac{1}{s_{2}^{fin}}} \right)}\phi_{f}}}}\end{matrix}$

As a matter of interest, the principal equivalent fracturepermeabilities obtained from the aforementioned calculations arek₁=k_(eq) ^(max) and k₂=k_(eq) ^(Min) in the principal directions 1 and2. The conductivities of the fractures of the fine equivalent network,C_(f1) ^(fin) and C_(f2) ^(fin), in these two directions of flow arededuced by writing the conservation of the flows per unit area offractured medium:

-   -   C_(f1) ^(fin)=s₂ ^(fin)·k₁ and C_(f2) ^(fin)=s₁ ^(fin)·k₂

Finally, the matrix medium between fractures has a permeability k_(m)^(fin).

-   -   Replacement of this fine network (FIG. 4) by the so-called        “coarse” equivalent network (FIG. 6) comprises more spaced-out        fractures to increase the degree of simplification with a view        to later flow simulations. The geometric and flow properties of        this coarse network are as follows:        -   fracture spacings s₁ ^(gros) and s₂ ^(gros) such that:    -   s₁ ^(gros)=G·s₁ ^(fin)    -   s₂ ^(gros)=G·s₂ ^(fin),

where G is a (fracture spacing) magnification coefficient of the networkwhose value is left to the user's discretion, with however an upperlimit G_(max-zone) that should not be exceeded to guarantee sufficientconnectivity between simplified and non-simplified zones:

G<G_(max-zone) where G_(max-zone) is such that

Max(s₁ ^(gros),s₂ ^(gros))<(minimum lateral dimension of the sub-zone)/6

Thus:

$G_{\max - {zone}} = \frac{DLM}{6 \cdot {{Max}\left( {s_{1}^{fin},s_{2}^{fin}} \right)}}$

with:

-   -   DLM is a minimum lateral dimension of the given sub-zone;    -   s₁ ^(fin) and s₂ ^(fin) are fracture spacings in the so-called        Warren and Root representation.    -   fracture conductivities C_(f1) ^(gros) and C_(f2) ^(gros) have        values allow keeping the flows per unit area of fractured        medium, that is also the equivalent permeabilities, that is:        C_(f1) ^(gros)=s₂ ^(gros)·k₁ and C_(f2) ^(gros)=s₁ ^(gros)·k₂,        or,        knowing that C_(f1) ^(fin)=s₂ ^(fin)·k₁ and C_(f2) ^(fin)=s₁        ^(fin)·k₂,        C_(f1) ^(gros)=G·C_(f1) ^(fin) and C_(f2) ^(gros)=G·C_(f2)        ^(fin)    -   A fracture opening e^(gros) allows again keeping the fracture        porosity φ_(f) of the initial network which is equal to that of        the coarse equivalent network:

ϕ_(f) = e^(gros) ⋅ (1/s₁^(gros) + 1/s₂^(gros)) = e^(fin) ⋅ (1/s₁^(fin) + 1/s₂^(fin))${{hence}\text{:}\mspace{14mu} e^{gros}} = {{e^{fin}\frac{\frac{1}{s_{1}^{fin}} + \frac{1}{s_{2}^{fin}}}{\frac{1}{s_{1}^{gros}} + \frac{1}{s_{2}^{gros}}}} = {G \cdot e^{fin}}}$

-   -   a matrix permeability k_(m) ^(gros) keeps the value of the        matrix-fracture exchange parameter:

$\begin{matrix}{\lambda_{fin} = {r_{w}^{2}\frac{k_{m}^{fin}}{k_{f}^{fin}}\left( {\frac{\alpha}{s_{1}^{{fin}^{2}}} + \frac{\alpha}{s_{2}^{{fin}^{2}}}} \right)}} \\{= \lambda_{grossier}} \\{= {r_{w}^{2}\frac{k_{m}^{gros}}{k_{f}^{gros}}\left( {\frac{\alpha}{s_{1}^{{gros}^{2}}} + \frac{\alpha}{s_{2}^{{gros}^{2}}}} \right)}}\end{matrix}$where α is a constant and where the fracture equivalent permeabilitiesof the fine and coarse networks are equal and denoted by k_(f) ^(fin)and k_(f) ^(gros),

That is:

$k_{m}^{gros} = {{k_{m}^{fin}\frac{\frac{1}{S_{1}^{{fin}^{2}}} + \frac{1}{S_{2}^{{fin}^{2}}}}{\frac{1}{S_{1}^{{gros}^{2}}} + \frac{1}{S_{2}^{{gros}^{2}}}}} = {G^{2}k_{m}^{fin}}}$

Finally, once this fine network coarse network→equivalence operationachieved for each sub-zone (or “block” of ZP), the fractures of thesimplified networks obtained are extended outside the limits ofsub-zones (or “blocks”) in order to guarantee sufficient partial“covering” and therefore sufficient horizontal connectivity of thesimplified networks of neighboring sub-zones (or “blocks” of ZP).Therefore, by following a test-proven procedure, the fractures of thesimplified network can thus be extended by a length equal to 60% of themaximum spacing (s₁ ^(gros),s₂ ^(gros)) of the fractures of thisnetwork.

Such an image is represented in FIG. 6.

Calibration of the Flow Properties of the Fractures

The next stage is the calibration of the flow properties of thefractures (fracture conductivity and opening), locally around the wells.This requires a well test simulation. According to the invention, thiswell test simulation is performed on the simplified flow models (FIG.6).

This type of calibration is well known. The method described in FrenchPatent 2,787,219 can for example be used. The flow responses of somewells (transient or pseudo-permanent flow tests, interferences, flowrate measurement, etc.) are simulated on these models extracted from thegeological model giving a discrete (realistic) representation of thefractures supplying these wells. The simulation result is then comparedwith the real measurements performed in the wells. If the resultsdiffer, the statistical parameters (PSF) describing the fracturenetworks are modified, then the flow properties of the initial fracturesare redetermined and a new simulation is carried out. The operation isrepeated until the simulation results and the measurements agree.

The results of these simulations allow calibration (estimation) of thegeometry and the flow properties of the fractures, such as theconductivities of the fracture networks of the reservoir being studiedand the openings.

3—Simulation of the Fluid Flows (SIM) and Optimization of the ReservoirProduction Conditions (OPT)

At this stage, the reservoir engineer has all the data required toconstruct the flow model on the reservoir scale. In fact, fracturedreservoir simulations often adopt the “double-porosity” approachproposed for example by Warren J. E. et al. in “The Behavior ofNaturally Fractured Reservoirs”, SPE Journal (September 1963) at pages245-255, according to which any elementary volume (cell of the reservoirmodel) of the fractured reservoir is modelled in a form of a set ofidentical parallelepipedic blocks, referred to as equivalent blocks,delimited by an orthogonal system of continuous uniform fracturesoriented in the principal directions of flow. The fluid flow on thereservoir scale occurs through the fractures only, and fluid exchangestake place locally between the fractures and the matrix blocks. Thereservoir engineer can for example use the methods described in thefollowing documents, applied to the entire reservoir this time: FrenchPatent 2,757,947 corresponding to U.S. Pat. No. 6,023,656 and FrenchPatent 2,757,957 corresponding to U.S. Pat. No. 6,064,944, and EP Patent2,037,080. These methods allow calculation of the equivalent fracturepermeabilities and the equivalent block dimensions for each cell of thereservoir model.

The reservoir engineer chooses a production process, for example thewaterflooding recovery process, for which the optimum implementationscenario remains to be specified for the field being considered. Thedefinition of an optimum waterflooding scenario is for example thesetting of the number and the location (position and spacing) of theinjector and producer wells in order to best account for the impact ofthe fractures on the progression of the fluids within the reservoir.

According to the scenario being selected, that is the double-mediumrepresentation of the reservoir and to the formula relating the massand/or energy exchange flow to the matrix-fracture potential difference,it is then possible to simulate the expected hydrocarbon production bymeans of the flow simulator (software) referred to as double-mediumsimulator.

At any time t of the simulated production, from input data E(t) (fixedor simulated-time varying data) and from the formula relating exchangeflow (f) to potential difference (ΔΦ), the simulator solves all theequations specific to each cell and each one of the two grids of themodel (equations involving the matrix-fracture exchange formuladescribed above), and it thus delivers the values of solution to theunknowns S(t) (saturations, pressures, concentrations, temperature,etc.) at this time t. This solution provides knowledge of the amounts ofoil produced and of the state of the reservoir (pressure distribution,saturations, etc.) at the time being considered.

4—Optimized Reservoir Development (EXPLO)

Selecting various scenarios characterized, for example, by variousrespective sites for the injector and producer wells, and simulating thehydrocarbon production for each one according to stage 3, enablesselection of the scenario allowing the production of the fracturedreservoir being considered to be optimized according to thetechnico-economic selected criteria.

Then the reservoir is developed according to this scenario allowing thereservoir production to be optimized.

The invention claimed is:
 1. A method for optimizing the development ofa fluid reservoir traversed by a fracture network and by at least onewell, wherein a representation of the fluid reservoir is constructed,the reservoir is discretized into a set of cells and the fractures arecharacterized by statistical parameters from observations of thereservoir, comprising: a) determining from the statistical parameters anequivalent permeability tensor and an average opening for the fracturesfrom which an image representative of the fracture network delimitingporous blocks and fractures is constructed; b) determining from thetensor a direction of flow of the fluid around the well; c) definingaround the well a first elliptical boundary zone centered on the welland containing the well and at least a second elliptical boundary zonecentered on the well which forms an elliptical ring with the ellipticalboundary of the first zone with the at least a second elliptical zonebeing oriented in the direction of flow of the fluid; d) simplifying theimage representative of the fracture network in a different manner ineach of the zones; e) using the simplified image to construct arepresentation of the fluid reservoir; and f) using the representationof the fluid reservoir and a flow simulator programmed with softwareexecuted by a computer to optimize the development of the fluidreservoir.
 2. A method as claimed in claim 1, wherein the statisticalparameters are selected from among the parameters: fracture density,fracture length, fracture orientation, fracture inclination, fractureopening and fracture distribution within the reservoir.
 3. A method asclaimed in claim 2, wherein an aspect ratio is determined for each zonewhich is defined from lengths of axes of an ellipse making up theboundary of each zone to reproduce a flow anisotropy around the wellwith the zones being constructed in accordance with the aspect ratio. 4.A method as claimed in claim 2, wherein a distance is defined betweenboundaries between the zones to give equal weight to each zone in termsof a pressure difference recorded in each zone under permanent flowconditions.
 5. A method as claimed in claim 2, comprising: repeating a)while modifying the statistical parameters to minimize a differencebetween a well test result and a well test simulation result from thesimplified image; associating with each one of the cells at least oneequivalent permeability value and an average opening value for thefractures with the values being determined from the modified statisticalparameters; simulating fluid flows in the reservoir by a flow simulatorequivalent permeability values and average opening values of thefractures associated with each of the cells; selecting a productionscenario optimizing the reservoir production by the fluid flowsimulation; and developing the reservoir according to the scenario tooptimize reservoir production.
 6. A method as claimed in claim 3,wherein the aspect ratio is determined by values of the permeabilitytensor.
 7. A method as claimed in claim 3, wherein a distance is definedbetween boundaries between the zones to give equal weight to each zonein terms of a pressure difference recorded in each zone under permanentflow conditions.
 8. A method as claimed in claim 4, wherein the distanceis defined by setting lengths of one of two axes of two successiveellipsis at values in geometric progression of a constant ratio.
 9. Amethod as claimed in claim 6, wherein a distance is defined betweenboundaries between the zones to give equal weight to each zone in termsof a pressure difference recorded in each zone under permanent flowconditions.
 10. A method as claimed in claim 7, wherein the distance isdefined by setting lengths of one of two axes of two successive ellipsisat values in geometric progression of a constant ratio.
 11. A method asclaimed in claim 9, wherein the distance is defined by setting lengthsof one of two axes of two successive ellipsis at values in geometricprogression of a constant ratio.
 12. A method as claimed in claim 1,wherein an aspect ratio is determined for each zone which is definedfrom lengths of axes of an ellipse making up the boundary of each zoneto reproduce a flow anisotropy around the well with the zones beingconstructed in accordance with the aspect ratio.
 13. A method as claimedin claim 12, wherein the aspect ratio is determined by values of thepermeability tensor.
 14. A method as claimed in claim 12, wherein adistance is defined between boundaries between the zones to give equalweight to each zone in terms of a pressure difference recorded in eachzone under permanent flow conditions.
 15. A method as claimed in claim12, comprising: repeating a) while modifying the statistical parametersto minimize a difference between a well test result and a well testsimulation result from the simplified image; associating with each oneof the cells at least one equivalent permeability value and an averageopening value for the fractures with the values being determined fromthe modified statistical parameters; simulating fluid flows in thereservoir by a flow simulator equivalent permeability values and averageopening values of the fractures associated with each of the cells;selecting a production scenario optimizing the reservoir production bythe fluid flow simulation; and developing the reservoir according to thescenario to optimize reservoir production.
 16. A method as claimed inclaim 13, wherein a distance is defined between boundaries between thezones to give equal weight to each zone in terms of a pressuredifference recorded in each zone under permanent flow conditions.
 17. Amethod as claimed in claim 13, comprising: repeating a) while modifyingthe statistical parameters to minimize a difference between a well testresult and a well test simulation result from the simplified image;associating with each one of the cells at least one equivalentpermeability value and an average opening value for the fractures withthe values being determined from the modified statistical parameters;simulating fluid flows in the reservoir by a flow simulator equivalentpermeability values and average opening values of the fracturesassociated with each of the cells; selecting a production scenariooptimizing the reservoir production by the fluid flow simulation; anddeveloping the reservoir according to the scenario to optimize reservoirproduction.
 18. A method as claimed in claim 14, wherein the distance isdefined by setting lengths of one of two axes of two successive ellipsisat values in geometric progression of a constant ratio.
 19. A method asclaimed in claim 16, wherein the distance is defined by setting lengthsof one of two axes of two successive ellipsis at values in geometricprogression of a constant ratio.
 20. A method as claimed in claim 1,wherein a distance is defined between boundaries between the zones togive equal weight to each zone in terms of a pressure differencerecorded in each zone under permanent flow conditions.
 21. A method asclaimed in claim 20, wherein the distance is defined by setting lengthsof one of two axes of two successive ellipsis values in geometricprogression of a constant ratio.
 22. A method as claimed in claim 20,comprising: repeating a) while modifying the statistical parameters tominimize a difference between a well test result and a well testsimulation result from the simplified image; associating with each oneof the cells at least one equivalent permeability value and an averageopening value for the fractures with the values being determined fromthe modified statistical parameters; simulating fluid flows in thereservoir by a flow simulator equivalent permeability values and averageopening values of the fractures associated with each of the cells;selecting a production scenario optimizing the reservoir production bythe fluid flow simulation; and developing the reservoir according to thescenario to optimize reservoir production.
 23. A method as claimed inclaim 21, comprising: repeating a) while modifying the statisticalparameters to minimize a difference between a well test result and awell test simulation result from the simplified image; associating witheach one of the cells at least one equivalent permeability value and anaverage opening value for the fractures with the values being determinedfrom the modified statistical parameters; simulating fluid flows in thereservoir by a flow simulator equivalent permeability values and averageopening values of the fractures associated with each of the cells;selecting a production scenario optimizing the reservoir production bythe fluid flow simulation; and developing the reservoir according to thescenario to optimize reservoir production.
 24. A method as claimed inclaim 1, wherein three zones are constructed with a first zonecontaining the well with no simplification of the image being provided,a second zone is constructed in contact with the first zone wherein afirst simplification of the image is carried out and a third zone isconstructed in contact with the second zone wherein a secondsimplification of the image is carried out with the secondsimplification being more significant than the first simplification. 25.A method as claimed in claim 24, wherein the second and third zones aredivided into sub-zones by the steps: dividing the second zone into anumber of sub-zones equal in number to a number of blocks of cellspresent in the second zone with a block of cells designating a verticalstack of cells; and dividing the third zone by the steps of dividingevery degree a boundary of the third zone with each degree defining 360arcs, defining a sub-zone by connecting end points of each of the arcsto a center of an ellipse forming the boundary, for each of thesub-zones, calculating an equivalent fracture permeability tensor fromwhich an orientation of the flows in the sub-zone is determined,comparing equivalent fracture permeability values and flow orientationbetween neighboring sub-zones, and grouping neighboring sub-zonestogether into a single sub-zone when a difference between thepermeability values is below a first threshold and when a differencebetween the flow orientations is below a second threshold.
 26. A methodas claimed in claim 24, comprising: repeating a) while modifying thestatistical parameters to minimize a difference between a well testresult and a well test simulation result from the simplified image;associating with each one of the cells at least one equivalentpermeability value and an average opening value for the fractures withthe values being determined from the modified statistical parameters;simulating fluid flows in the reservoir by a flow simulator equivalentpermeability values and average opening values of the fracturesassociated with each of the cells; selecting a production scenariooptimizing the reservoir production by the fluid flow simulation; anddeveloping the reservoir according to the scenario to optimize reservoirproduction.
 27. A method as claimed in claim 25 wherein, for a givensub-zone, value Gmax-zone is equal to:$G_{\max - {zone}} = \frac{DLM}{6 \cdot {{Max}\left( {s_{1}^{fin},s_{2}^{fin}} \right)}}$with: DLM being a minimum lateral dimension of the given sub-zone; ands₁ ^(fin), s₂ ^(fin) being fracture spacings in the Warren and Rootrepresentation.
 28. A method as claimed in claim 25, comprising:repeating a) while modifying the statistical parameters to minimize adifference between a well test result and a well test simulation resultfrom the simplified image; associating with each one of the cells atleast one equivalent permeability value and an average opening value forthe fractures with the values being determined from the modifiedstatistical parameters; simulating fluid flows in the reservoir by aflow simulator equivalent permeability values and average opening valuesof the fractures associated with each of the cells; selecting aproduction scenario optimizing the reservoir production by the fluidflow simulation; and developing the reservoir according to the scenarioto optimize reservoir production.
 29. A method as claimed in claim 1,wherein the image is simplified by the steps: constructing a fracturenetwork equivalent to the image, by a Warren and Root representationwherein the network has fracture spacings in two orthogonal directionsof principal permeability, by a fracture opening parameter, by fractureconductivities and a permeability kmfin of a matrix medium betweenfractures; simplifying the equivalent fracture network by a networkfracture spacing coefficient whose value is less than a value Gmax-zonedefined for each of the zones to guarantee connectivity betweensimplified zones and non-simplified zones.
 30. A method as claimed inclaim 29, comprising: repeating a) while modifying the statisticalparameters to minimize a difference between a well test result and awell test simulation result from the simplified image; associating witheach one of the cells at least one equivalent permeability value and anaverage opening value for the fractures with the values being determinedfrom the modified statistical parameters; simulating fluid flows in thereservoir by a flow simulator equivalent permeability values and averageopening values of the fractures associated with each of the cells;selecting a production scenario optimizing the reservoir production bythe fluid flow simulation; and developing the reservoir according to thescenario to optimize reservoir production.
 31. A method as claimed inclaim 1, comprising: repeating a) while modifying the statisticalparameters to minimize a difference between a well test result and awell test simulation result from the simplified image; associating witheach one of the cells at least one equivalent permeability value and anaverage opening value for the fractures with the values being determinedfrom the modified statistical parameters; simulating fluid flows in thereservoir by a flow simulator equivalent permeability values and averageopening values of the fractures associated with each of the cells;selecting a production scenario optimizing the reservoir production bythe fluid flow simulation; and developing the reservoir according to thescenario to optimize reservoir production.